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Binary black-hole systems are expected to be important sources of gravitational waves for upcom¬ 
ing gravitational-wave detectors. If the spins are not colinear with each other or with the orbital 
angular momentum, these systems exhibit complicated precession dynamics that are imprinted on 
the gravitational waveform. We develop a new procedure to match the precession dynamics com¬ 
puted by post-Newtonian (PN) theory to those of numerical binary black-hole simulations in full 
general relativity. For numerical relativity (NR) simulations lasting approximately two precession 
cycles, we find that the PN and NR predictions for the directions of the orbital angular momentum 
and the spins agree to better than ~ 1° with NR during the inspiral, increasing to 5° near merger. 

Nutation of the orbital plane on the orbital time-scale agrees well between NR and PN, whereas 
nutation of the spin direction shows qualitatively different behavior in PN and NR. We also examine 
how the PN equations for precession and orbital-phase evolution converge with PN order, and we 
quantify the impact of various choices for handling partially known PN terms. 


I. INTRODUCTION 

Binary black holes (BBH) are among the most im¬ 
portant sources of gravitational waves for upcoming 
gravitational-wave detectors like Advanced LIGO [1] and 
Virgo [2]. Accurate predictions of the gravitational wave¬ 
forms emitted by such systems are important for detec¬ 
tion of gravitational waves and for parameter estimation 
of any detected binary [3] . When either black hole carries 
spin that is not aligned with the orbital angular momen¬ 
tum, there is an exchange of angular momentum between 
the components of the system, leading to complicated 
dynamical behavior. Figure 1 exhibits the directions of 
the various angular momenta in several simulations de¬ 
scribed in this paper. This behavior is imprinted on the 
emitted waveforms [4-6] , making them more feature-rich 
than waveforms from aligned-spin BBH systems or non¬ 
spinning BBH systems. In order to model the waveforms 
accurately, then, we need to understand the dynamics. 

The orbital-phase evolution of an inspiraling binary, 
the precession of the orbital angular momentum and 
the black-hole spins, and the emitted gravitational wave¬ 
forms can be modeled with post-Newtonian theory [7], 
a perturbative solution of Einstein’s equations in pow¬ 
ers of v/c, the ratio of the velocity of the black holes to 
the speed of light. Such post-Newtonian waveforms play 
an important role in the waveform modeling for ground- 
based interferometric gravitational-wave detectors (see, 
e.g., [8]). For non-spinning and aligned-spin BBH, how¬ 
ever, the loss of accuracy of the post-Newtonian phase 
evolution in the late inspiral has been identified as one of 
the dominant limitations of waveform modeling [9-14]. 

Precessing waveform models (e.g., [6, 15-17]) depend 
on the orbital phase evolution and the precession dynam¬ 
ics. Therefore, it is important to quantify the accuracy of 
the post-Newtonian approximation for modeling the pre¬ 


cession dynamics itself, and the orbital-phase evolution of 
precessing binaries. Recently, the SXS collaboration has 
published numerical-relativity solutions to the full Ein¬ 
stein equations for precessing BBH systems [18]. These 
simulations cover > 30 orbits and up to two precession 
cycles. Therefore, they offer a novel opportunity to sys¬ 
tematically quantify the accuracy of the post-Newtonian 
precession equations, the topic of this paper. 

In this paper, we develop a new technique to match 
the initial conditions of post-Newtonian dynamics to a 
numerical relativity simulation. We then use this tech¬ 
nique to study the level of agreement between the post- 
Newtonian precession equations and the numerical sim¬ 
ulations. The agreement is remarkably good, the direc¬ 
tions of orbital angular momentum and spin axes in post- 
Newtonian theory reproduces the numerical simulations 
usually to better than 1 degree. We also investigate nu¬ 
tation effects on the orbital time-scale that are imprinted 
both in the orbital angular momentum and the spin- 
directions. For the orbital angular momentum, NR and 
PN yield very similar nutation features, whereas for the 
spin direction, nutation is qualitatively different in PN 
and the investigated NR simulations. Considering the 
orbital-phase evolution, we find that the disagreement 
between post-Newtonian orbital phase and numerical rel¬ 
ativity simulation is comparable to the aligned-spin case. 
This implies that the orbital phase evolution will remain 
an important limitation for post-Newtonian waveforms 
even in the precessing case. Finally, we study the conver¬ 
gence with post-Newtonian order of the precession equa¬ 
tions, and establish very regular and fast convergence, in 
contrast to post-Newtonian orbital phasing. 

This paper is organized as follows: Section II describes 
the post-Newtonian expressions utilized, the numerical 
simulations, how we compare PN and NR systems with 
each other, and how we determine suitable “best-fitting” 
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FIG. 1. Precession cones of the six primary precessing simulations considered here, as computed by NR and PN. Shown are 
the paths traced on the unit sphere by the normal to the orbital plane i and the spin-directions Xi, 2 - The thick lines represent 
the NR data, with the filled circles indicating the start of the NR simulations. The lines connecting the NR data to the origin 
are drawn to help visualize the precession-cones. The PN data, plotted with thin lines, lie on the scale of this figure almost 
precisely on top of the NR data. (The PN data was constructed using the Taylor T4 approximant matched at frequency 
m£l m = 0.021067, with a matching interval width = 0.1f2 m .) 


PN parameters for a comparison with a given NR sim¬ 
ulation. Section III presents our results, starting with 
a comparison of the precession dynamics in Sec. Ill A, 
and continuing with an investigation in the accuracy of 
the orbital phasing in Sec. IIIB. The following two sec¬ 
tions study the convergence of the PN precession equa¬ 
tions and the impact of ambiguous choices when dealing 
with incompletely known spin-terms in the PN orbital 
phasing. Section HIE, finally, is devoted to some tech¬ 
nical numerical aspects, including an investigation into 
the importance of the gauge conditions used for the NR 
runs. We close with a discussion in Sec. IV. The appen¬ 
dices collect the precise post-Newtonian expressions we 
use and additional useful formulae about quaternions. 

II. METHODOLOGY 
A. Post-Newtonian Theory 

Post-Newtonian (PN) theory is an approximation to 
General Relativity in the weak-field, slow-motion regime, 
characterized by the small parameter e ~ (u/c) 2 ~ 


where in, v, and r denote the characteristic mass, veloc¬ 
ity, and size of the source, c is the speed of light, and G 
is Newton’s gravitational constant. For the rest of this 
paper, the source is always a binary black-hole system 
with total mass m, relative velocity v and separation r, 
and we use units where G = c = 1. 

Restricting attention to quasi-spherical binaries in the 
adiabatic limit, the local dynamics of the source can be 
split into two parts: the evolution of the orbital fre¬ 
quency, and the precession of the orbital plane and the 
spins. The leading-order precessional effects [19] and 
spin contributions to the evolution of the orbital fre¬ 
quency [20, 21] enter post-Newtonian dynamics at the 
1.5 PN order (i.e., e 3 / 2 ) for spin-orbit effects, and 2 PN 
order for spin-spin effects. We also include non-spin 
terms to 3.5 PN order [7], the spin-orbit terms to 4 PN 
order [22], spin-spin terms to 2 PN order [21] 1 . For the 
precession equations, we include the spin-orbit contribu- 


1 During the preparation of this manuscript, the 3 PN spin-spin 

contributions to the flux and binding energy were completed in 

[23]. These terms are not used in the analysis presented here. 
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FIG. 2. Vectors describing the orbital dynamics of the system. 
The yellow plane denotes the orbital plane. Rt(t) is the rotor 
that rotates the coordinate triad (x, y, z) into the orbital triad 

(ft, Ac¬ 


tions to next-to-next-to-leading order, corresponding to 
3.5 PN [24], The spin-spin terms are included at 2 PN 
order 2 . 


1. Orbital dynamics 

Following earlier work (e.g., Ref. [21]) we describe the 
processing BH binary by the evolution of the orthonormal 
triad (ft, A,£), as indicated in Fig. 2: ft denotes the unit 
separation vector between the two compact objects, t is 
the normal to the orbital plane and A = l x ft completes 
the triad. This triad is time-dependent, and is related to 
the constant inertial triad {x, y, z) by a time-dependent 
rotation Rf, as indicated in Fig. 2. The rotation Rf 
will play an important role in Sec. IIC. The orbital triad 
obeys the following equations: 


d[ 

dt 

dfi 

dt 

dX 

dt 


vuh x l, 

(la) 

fix, 

(lb) 

—flh + mt. 

(lc) 


2 The investigation of the effects of spin-spin terms at higher 
PN orders (see e.g. [25-27] and references therein), and terms 
which are higher order in spin (e.g cubic spin terms) [28, 29] is 
left for future work. 


Here, f l is the instantaneous orbital frequency and vj is 
the precession frequency of the orbital plane. 

The dimensionless spin vectors \i = Si/ntf also obey 
precession equations: 


d ^=Q 1 xxi, (2a) 

at 

^ = ^2 X X2- (2b) 

dt 


The precession frequencies fir 2 j vo are series in the PN 
expansion parameter e; their explicit form is given in Ap¬ 
pendix A. 

The evolution of the orbital frequency is derived from 
energy balance: 


dE 

dt 




( 3 ) 


where E is the energy of the binary and T is the 
gravitational-wave flux. E and T are PN series depend¬ 
ing on the orbital frequency fl, the vector £, and the 
BH spins X\ , X'l- Their explicit formulas are given in 
Appendix A. In terms of x = (mfl) 2 / 3 ~ e, Eq. (3) be¬ 
comes: 


dx T 

dt dE/dx ’ 


( 4 ) 


where the right-hand side is a ratio of two PN series. 

There are several well known ways of solving Eq. (4), 
which lead to different treatment of uncontrolled higher- 
order PN terms—referred to as the Taylor T1 through T5 
approximants [30, 31]. The T2 and T3 approximants can¬ 
not be applied to general precessing systems; we there¬ 
fore exclude them from this work. We now briefly review 
the remaining approximants, which will be used through¬ 
out this work. ’ The most straightforward approach is 
to evaluate the numerator and denominator of Eq. (4) 
and then solve the resulting ordinary differential equa¬ 
tion numerically, which is the Taylor T1 approximant. 
Another approach is to re-expand the ratio F / {dE / dx) 
in a new power series in x, and then truncate at the ap¬ 
propriate order. This gives the Taylor T4 approximant. 
Finally, one can expand the inverse of the right-hand-side 
of Eq. (4) in a new power series in x, truncate it at the 
appropriate order, and then substitute the inverse of the 
truncated series into the right-hand side in Eq. (4) . This 
last approach, known as the Taylor T5 approximant [31], 
has been introduced fairly recently. 


2. Handling of spin terms 

When constructing Taylor approximants that include 
the re-expansion of the energy balance equation, the han¬ 
dling of spin terms becomes important. In particular, 


3 See, e.g., Ref. [32] for a more complete description of approxi¬ 
mants T1 through T4. 
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TABLE I. Numerical relativity simulations utilized here. SXS ID refers to the simulation number in Ref. [18], q = m\/m 2 
is the mass ratio, yi .2 are the dimensionless spins, given in coordinates where h(t = O') = x, £(t = 0) = z. Do, flo and e are 
the initial coordinate separation, the initial orbital frequency, and the orbital eccentricity, respectively. The first block lists 
the precessing runs utilized, where \i,r = (—0.18,-0.0479,-0.0378) and \ 2 ,r = (-0.0675,0.0779,-0.357). The second block 
indicates 31 further precessing simulations used in Fig. 10, and the last block lists the aligned spin systems for orbital phase 
comparisons. 


Name 

SXS ID 

q 

Xi 

X2 

Do/M 

?nflo 

e 

ql_0.5x 

0003 

1.0 

(0.5,0.0,0) 

(0,0,0) 

19 

0.01128 

0.003 

ql.5_0.5x 

0017 

1.5 

(0.5,0,0) 

(0,0,0) 

16 

0.01443 

< 2 x 10 -4 

q3_0.5x 

0034 

3.0 

(0.5,0,0) 

(0,0,0) 

14 

0.01743 

< 2 x 10“ 4 

q5_0.5x 


5.0 

(0.5,0,0) 

(0,0,0) 

15 

0.01579 

0.002 

ql_two_spins 

0163 

1.0 

(0.52,0,-0.3) 

(0.52,0,0.3) 

15.3 

0.01510 

0.003 

ql.97_random 

0146 

1.97 

Xl ,r 

X2,r 

15 

0.01585 

1 

O 

V 

31 random runs 

115-145 

[R2] 

Xi < 0.5 

X2 < 0.5 

15 

« 0.0159 

[10- 4 ,10~ 3 ] 

ql_0.5z 

0005 

1.0 

(0,0,0.5) 

(0,0,0) 

19 

0.01217 

0.0003 

ql_-0.5z 

0004 

1.0 

(0,0,0.5) 

(0,0,0) 

19 

0.01131 

0.0004 

ql.5_0.5z 

0013 

1.5 

(0,0,0.5) 

(0,0,0) 

16 

0.01438 

0.00014 

ql.5_-0.5z 

0012 

1.5 

(0,0,-0.5) 

(0,0,0) 

16 

0.01449 

0.00007 

q3_0.5z 

0031 

3.0 

(0,0,0.5) 

(0,0,0) 

14 

0.01734 

< 10” 4 

q3_-0.5z 

0038 

3.0 

(0,0,-0.5) 

(0,0,0) 

14 

0.01756 

< 10" 4 

q5_0.5z 

0061 

5.0 

(0,0,0.5) 

(0,0,0) 

15 

0.01570 

0.004 

q5_-0.5z 

0060 

5.0 

(0,0,-0.5) 

(0,0,0) 

15 

0.01591 

0.003 

q8_0.5z 

0065 

8.0 

(0,0,0.5) 

(0,0,0) 

13 

0.01922 

0.004 

q8_-0.5z 

0064 

8.0 

(0,0,-0.5) 

(0,0,0) 

13 

0.01954 

0.0005 


terms of quadratic and higher order in spins, such as 
(Si) 2 , appear in the evolution of the orbital frequency 
at 3 PN and higher orders. These terms arise from 
lower-order effects and represent incomplete information, 
since the corresponding terms are unknown in the orig¬ 
inal power series for the binding energy E and the flux 

r, 


E(x) 

F(x) 



( 5 ) 

( 6 ) 


where m — m\ + m 2 and v = rtixmi/rn 2 , and 7714,2 are 
the individual masses. 

In these expansions, the spin-squared terms come in at 
2 PN order and thus appear in 04 and 64 , cf. Eqs. (A18) 
and (A24). Then, in the re-expansion series of Taylor 
T4, 


S 


T 

dE/dx 


64ZZ r 

5m ^ 


(1 + ^2s k x k/2 ), 

k—2 


( 7 ) 


the coefficients Sk can be recursively determined, e.g. 


S4 = &4 — 3G4 — 2S2<22, (8) 

5 

se = b 6 - (4a 6 + 3 s 2 a 4 + -s 3 a 3 + 2s 4 a 2 ). (9) 

Thus, the spin-squared terms in 04 and 64 will induce 
spin-squared terms at 3PN order in S 6 - The analogous 


conclusion holds for Taylor T5. These spin-squared terms 
are incomplete as the corresponding terms in the binding 
energy and flux (i.e. in ae and be) are not known. 

This re-expansion has been handled in several ways 
in the literature. For example, Nitz et al. [14] include 
only terms which are linear in spin beyond 2 PN order. 
On the other hand, Santamaria et al. [33] keep all terms 
in spin arising from known terms in E and T. In the 
present work, we also keep all terms up to 3.5 PN order, 
which is the highest order to which non-spin terms are 
completely known. Similarly, we include all terms when 
computing the precession frequency (see A 2) . We inves¬ 
tigate the impact of different spin-truncation choices in 
Sec. Ill D, along with the impact of partially known 4 PN 
spin terms. 


B. Numerical Relativity Simulations 

To characterize the effectiveness of PN theory in re¬ 
producing NR results, we have selected a subset of 16 
simulations from the SXS waveform catalog described in 
Ref. [18]. 4 Our primary results are based on six precess¬ 
ing simulations and a further ten non-precessing ones 
for cross-comparisons. To check for systematic effects, 
we use a further 31 precessing simulations with random 
mass-ratios and spins. The parameters of these runs are 


4 The waveform and orbital data are publicly available at https: 
//www.black-holes.org/waveforms/. 
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given in Table I. They were chosen to represent vari¬ 
ous degrees of complexity in the dynamics: (i) precessing 
versus non-precessing simulations, the latter with spins 
parallel or anti-parallel to t\ (ii) one versus two spinning 
black holes; (iii) coverage of mass ratio from q = 1 to 
q = 8; (iv) long simulations that cover more than a pre¬ 
cession cycle; and (v) a variety of orientations of yq, X 2 , £■ 
Figure 1 shows the precession cones of the normal to the 
orbital plane and the spins for for the six primary pre¬ 
cessing cases in Table I. The PN data were computed 
using the Taylor T4 3.5 PN approximant. 

The simulations from the catalog listed in Table I were 
run with numerical methods similar to [34]. A gener¬ 
alized harmonic evolution system [35-38] is employed, 
and the gauge is determined by gauge source functions 
H a . During the inspiral phase of the simulations consid¬ 
ered here, H a is kept constant in the co-moving frame, 
cf. [32, 39, 40]. About 1.5 orbits before merger, the gauge 
is changed to damped harmonic gauge [41-43]. This 
gauge change happens outside the focus of the compar¬ 
isons presented here. 

The simulation q5_0.5x analyzed here is a re-run of the 
SXS simulation SXS:BBH:0058 from Ref. [18]. We per¬ 
formed this re-run for two reasons: First, SXS:BBH:0058 
changes to damped harmonic gauge in the middle of the 
inspiral, rather than close to merger as all other cases 
considered in this work. Second, SXS:BBH:0058 uses an 
unsatisfactorily low numerical resolution during the cal¬ 
culation of the black hole spins. Both these choices leave 
noticeable imprints on the data from SXS:BBH:0058, and 
the re-run q5_0.5x allows us to quantify the impact of 
these deficiencies. We discuss these effects in detail in 
Secs. HIE2 and HIE3. The re-run q5_0.5x analyzed 
here is performed with improved numerical techniques. 
Most importantly, damped harmonic gauge is used essen¬ 
tially from the start of the simulation, t > 100 M. The 
simulation q5_0.5x also benefits from improved adaptive 
mesh refinement [44] and improved methods for control¬ 
ling the shape and size of the excision boundaries; the 
latter methods are described in Sec.II.B. of Ref. [45]. 

We have performed convergence tests for some of the 
simulations; Sec. HIE will demonstrate with Fig. 16 that 
numerical truncation error is unimportant for the com¬ 
parisons presented here. 


C. Characterizing Precession 

The symmetries of non-precessing systems greatly sim¬ 
plify the problem of understanding the motion of the bi¬ 
nary. In a non-precessing system, the spin vectors are 
essentially constant, and two of the rotational degrees of 
freedom are eliminated in the binary’s orbital elements. 
Assuming quasi-circular orbits, the entire system can be 
described by the orbital phase <E>, which can be defined 
as the angle between h and x. In post-Newtonian theory 
the separation between the black holes can be derived 
from d^/dt. Thus comparison between post-Newtonian 


and numerical orbits, for example, reduces entirely to 
the comparison between <I>pn and 4 >nr [32, 46]. For pre¬ 
cessing systems, on the other hand, the concept of an 
orbital phase is insufficient; *I> could be thought of as 
just one of the three Euler angles. We saw in Sec. IIA 1 
that the orbital dynamics of a precessing system can be 
fairly complex, involving the triad (n, A, t) (or equiva¬ 
lently the frame rotor Rf) as well as the two spin vectors 
Xi and \ 2 —each of which is, of course, time dependent. 
When comparing post-Newtonian and numerical results, 
we need to measure differences between each of these 
quantities in their respective systems. 

To compare the positions and velocities of the black 
holes themselves, we can condense the information about 
the triads into the quaternion quantity [47] 

Ra ■= Rf N Rf R , (io) 

which represents the rotation needed to align the PN 
frame with the NR frame. This is a geometrically mean¬ 
ingful measure of the relative difference between two 
frames. We can reduce this to a single real number by 
taking the magnitude of the logarithm of this quantity, 
defining the angle 5 

:=2|logJ? A | ■ (11) 

This measure has various useful qualities. It is invariant, 
in the sense that any basis frame used to define I? PN and 
Rf 111 will result in the same value of $a- It conveniently 
distills the information about the difference between the 
frames into a single value, but is also non-degenerate in 
the sense that $a = 0 if and only if the frames are iden¬ 
tical. It also reduces precisely to <f> PN — $ NR for non- 
precessing systems; for precessing systems it also incor¬ 
porates contributions from the relative orientations of the 
orbital planes. 6 

Despite these useful features of $a, it may sometimes 
be interesting to use different measures, to extract indi¬ 
vidual components of the binary evolution. For example, 
Eq. (la) describes the precession of the orbital plane. 
When comparing this precession for two approaches, a 
more informative quantity than 4 >a is simply the angle 
between the £ vectors in the two systems: 

ZL = cos" 1 (A N • ^ NR ) . (12) 

Similarly, we will be interested in understanding the evo¬ 
lution of the spin vectors, as given in Eqs. (2). For this 
purpose, we define the angles between the spin vectors: 

^Xi — cos -1 (x™ ' Xi* R ) i (13a) 


5 More explanation of these expressions, along with relevant for¬ 
mulas for calculating their values, can be found in Appendix B. 

6 It is interesting to note that any attempt to define the orbital 

phases of precessing systems separately, and then compare them 
as some is either ill defined or degenerate—as shown 

in Appendix B 4. This does not mean that it is impossible to 
define such phases, but at best they will be degenerate; multiple 
angles would be needed to represent the full dynamics. 



6 



t - imerge (M) 



t ~ tmerge (M) 

FIG. 3. Examples of the averaging procedure and error 
estimates employed for all comparisons. Shown here are 
ql.97_random and q5.0_0.5x. PN evolutions were performed 
with the Taylor T1 approximant. The thin blue lines show all 
the PN-NR matching intervals. 

ZX 2 = cos- 1 • xD • (13b) 

We will use all four of these angles below to compare the 
post-Newtonian and numerical orbital elements. 

D. Matching Post-Newtonian to Numerical 
Relativity 

When comparing PN theory to NR results, it is im¬ 
portant to ensure that the initial conditions used in both 
cases represent the same physical situation. We choose a 
particular orbital frequency fi m and use the NR data to 
convert it to a time t rn . To initialize a PN evolution at 
t rn . we need to specify 

Q,Xi,X2, (14) 

t,n,X.uX 2 , (15) 

fi. (16) 

The quantities (14) are conserved during the PN evolu¬ 
tion. The quantities (15) determine the orientation of 


the the binary and its spins relative to the inertial triad 
(x,y,z). The orbital frequency 14 in Eq. (16), finally, 
parametrizes the separation of the binary at t m . The 
simplest approach is to initialize the PN evolution from 
the respective quantities in the initial data of the NR evo¬ 
lution. This would neglect initial transients in NR data 
as in, e.g., Fig. 1 of Ref. [40]. These transients affect the 
masses and spins of the black holes, so any further PN- 
NR comparisons would be comparing slightly different 
physical configurations. The NR transients decay away 
within the first orbit of the NR simulation, so one can 
consider initializing the PN evolution from NR at a time 
after the NR run has settled down. However, the gen¬ 
erally non-zero (albeit very small) orbital eccentricity in 
the NR simulation can lead to systematic errors in the 
subsequent comparison as pointed out in Ref. [32]. 

Therefore, we use time-averaged quantities evaluated 
after the initial transients have vanished. In particular, 
given a numerical relativity simulation, we set the PN 
variables listed in Eq. (14) to their numerical relativity 
values after junk radiation has propagated away. 

The remaining nine quantities Eqs. (15) and (16) must 
satisfy the constraint i-n = 0. We determine them with 
constrained minimization by first choosing an orbital fre¬ 
quency interval [fi m — 6CI/2, fi m + <5fi/2] of width Sfl. 
Computing the corresponding time interval [ti,tf] in the 
NR simulation, we define the time average of any quan¬ 
tity Q by 

(Q) = r^-r I ' 1 Qdt. (17) 

tf ~ H J ti 

Using these averages, we construct the objective func¬ 
tional S as 

5 = {{XL) 2 ) + <(Z X i) 2 > + <(Z X2 ) 2 ) + ((AH) 2 ) (18) 

where AH = (fip^ — Hnr)/Hnr- When a spin on 
the black holes is below 10” 5 the corresponding term is 
dropped from Eq. (18). The objective functional is then 
minimized using the SLSQP algorithm [48, 49] to allow 
for constrained minimization. In Eq. (18) we use equal 
weights for each term; other choices of the weights do not 
change the qualitative picture that we present. 

The frequency interval [fi m ± Jfi/2] is chosen based 
on several considerations. First it is selected after junk 
radiation has propagated away. Secondly, it is made wide 
enough so that any residual eccentricity effects average 
out. Finally, we would like to match PN and NR as early 
as possible. But since we want to compare various cases 
to each other, the lowest possible matching frequency will 
be limited by the shortest NR run (case q8_-0.5z). Within 
these constraints, we choose several matching intervals, 
in order to estimate the impact of the choice of matching 
interval on our eventual results. Specifically, we use three 
matching frequencies 

rnfi m e {0.021067,0.021264,0.021461}, (19) 
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and employ four different matching windows for each, 
namely 

sn/n rn e {0.06,0.08,o.i, 0.12}. (20) 

These frequencies correspond approximately to the range 
between 10-27 orbits to merger depending on the param¬ 
eters of the binary, with the lower limit for the case ql.0_- 
0.5x and the upper for q8.0_0.5x. 

Matching at multiple frequencies and frequency win¬ 
dows allows an estimate on the error in the matching 
and also ensures that the results are not sensitive to the 
matching interval being used. In this article, we gener¬ 
ally report results that are averaged over the 12 PN-NR 
comparisons performed with the different matching in¬ 
tervals. We report error bars ranging from the smallest 
to the largest result among the 12 matching intervals. 
As examples, Fig. 3 shows <E>a as a function of time to 
merger f m erge for the cases ql.97_random and q5_0.5x for 
all the matching frequencies and intervals, as well as the 
average result and an estimate of the error. Here t melge is 
the time in the NR simulation when the common horizon 
is detected. 



t Emerge (M) 


FIG. 4. Angle ZL by which f PN (t) differs from i NR (t) for 
the configuration ql_0.5x (red lines) and q5_0.5x (black lines). 
ZL < 0.2° except very close to merger. In each case, the PN 
predictions based on different PN approximants are shown 
in different line styles. Shown is the point-wise average of 
12 ZL(t) curves, i.e. the thick red line of Fig. 3. The thin 
horizontal lines show the widest edges of the PN matching 
intervals. 


III. RESULTS 
A. Precession Comparisons 

We apply the matching procedure of Sec. IID to the 
precessing NR simulations in Table I. PN-NR matching 
is always performed at the frequencies given by Eq. (19) 
which are the lowest feasible orbital frequencies across all 
cases in Table I. Figure 1 shows the precession cones for 
the normal to the orbital plane l and the spins Xi, 2 - As 
time progresses, t and yq .2 undergo precession and nuta¬ 
tion, and the precession cone widens due to the emission 
of gravitational radiation. Qualitatively, the PN results 
seem to follow the NR results well, until close to merger. 

We now turn to a quantitative analysis of the preces¬ 
sion dynamics, establishing first that the choice of Taylor 
approximant is of minor importance for the precession 
dynamics. We match PN dynamics to the NR simula¬ 
tions q5_0.5x and ql_0.5x for the Taylor approximants 
Tl, T4 and T5. We then compute the angles ZL and 
Zyq. Figure 4 shows the resulting ZL. During most of 
the inspiral, we find ZL of order a few 10 -3 radians in¬ 
creasing to ~ 0.1 radians during the last 1000M before 
merger. Thus the direction of the normal to the orbital 
plane is reproduced well by PN theory. This result is vir¬ 
tually independent of the Taylor approximant suggesting 
that the choice of approximant only weakly influences 
how well PN precession equations track the motion of 
the orbital plane. In other words, precession dynamics 
does not depend on details of orbital phasing like the 
unmodeled higher-order terms in which the Taylor ap¬ 
proximants differ from each other. 

Turning to the spin direction yq we compute the angle 
Zyq between y;^ R (f) an d Xi N (t) and plot the result in 
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FIG. 5. Angle Zyi by which xf N (t) differs from y{ R (t) for the 
configuration ql_0.5x (red lines) and q5_0.5x (black lines). In 
each case, the PN predictions based on different PN approxi¬ 
mants are shown in different line styles. The thin horizontal 
lines show the widest edges of the PN matching intervals. 

Fig. 5. While Fig. 5 looks busy, the first conclusion is 
that Zyi is quite small < 0.01 rad through most of the 
inspiral, and rises somewhat close to merger. 

The pronounced short-period oscillations of Zy 1 in 
Fig. 5 are caused by differences between PN-nutation fea¬ 
tures and NR-nutation features. To better understand 
the nutation features and their impact on the angle Zyq, 
we remove nutation features by filtering out all frequen¬ 
cies comparable to the orbital frequency. This is possible 
because the precession frequency is much smaller than 
the nutation frequency. The filtering is performed with a 
3rd order, bi-directional low pass Butterworth filter [50] 
with a fixed cutoff frequency chosen to be lower than the 









FIG. 6. Angle Z\ 1 between the “orbit-averaged” spins for 
the configuration q5_0.5x. The non orbit-averaged difference 
Z\i (cf. Fig. 5) is shown for comparison. Shown is one match¬ 
ing interval as indicated by the thin horizontal line. 

nutation frequency at the start of the inspiral. Due to 
the nature of the filtering, the resulting averaged spin 
will suffer from edge effects which affect approximately 
the first and last 1000 M of the inspiral. Furthermore, 
the precession frequency close to merger becomes compa¬ 
rable to the nutation frequency at the start of the simula¬ 
tion and thus filtering is no longer truthful in this region. 
Therefore, we only use the “averaged” spins where such 
features are absent. 

Applying this smoothing procedure to both xf N and 
^nr £ or ^e run q 5 _Q_ 5 Xj we compute the angle Z\i be¬ 
tween the averaged spin vectors, x PN and x^ R . This an¬ 
gle is plotted in Fig. 6 7 , where results only for the Taylor 
T1 approximant are shown, and for only one matching in¬ 
terval specified by mil rn = 0.0210597 and 8fl/n m = 0.1. 
The orbit-averaged spin directions ^N R / PN agree signifi¬ 
cantly better with each other than the non-averaged ones 
(cf. the black line in Fig. 6, which is duplicated from 
Fig. 5). In fact, the orbit-averaged spin precessing be¬ 
tween NR and PN agrees as well as the orbital angular 
momentum precession, cf. Fig. 4. Thus, the difference in 
the spin dynamics is dominated by the nutation features, 
with the orbit-averaged spin dynamics agreeing well be¬ 
tween PN and NR. 

To characterize the nutation features in the spin vec¬ 
tors, we introduce a coordinate system which is specially 
adapted to highlighting nutation effects. The idea is to 
visualize nutation with respect to the averaged spin vec¬ 
tor x■ We compute the time-derivative x numerically. 
Assuming that the “averaged” spin is undergoing pure 
precession, so that x • X = 0, we define a new coordinate 
system (ei,e 2 ,e 3 ) by ei = x, e 2 = x/|x|,e 3 = q x e 2 . 
The spin is now projected onto the e 2 — e 3 plane, thus 


' To illustrate edge effects of the Butterworth filter, Fig. 6 includes 
the early and late time periods where the filter affects Z \i ■ 



FIG. 7. The projection of Xi™ and xf N onto the e 2 — e 3 
plane described in the text for case q5_0.5x. The system is 
shown in the interval t — t me rge £ [—6662, —1556]. along the 
e 3 axis. Meanwhile, the NR data show variations in e 2 and 
e 3 directions of comparable magnitude. The solid symbols 
(black diamond for NR, red square for PN) indicate the data 
at the start of the plotted interval, chosen such that ' n 
is maximal—i.e., where the spin projection into the orbital 
plane is parallel to n. The subsequent four open symbols 
(blue diamonds for NR, green squares for PN) indicating the 
position 1/8-th, 1/4-th, 3/8-th and 1/2 of an orbit later. 


showing the motion of the spin in a frame “coprecessing” 
with the averaged spin. This allows us to approximately 
decouple precession and nutation and compare them sep¬ 
arately between PN and NR. 

Figure 7 plots the projection of the spins Xi™ an d xf N 
onto their respective “orbit averaged” e 2 — e 3 planes. We 
see that the behavior of the NR spin and the PN spins are 
qualitatively different: For this single-spin system, the 
PN spin essentially changes only in the e 3 direction (i.e., 
orthogonal to its average motion x PN ). In contrast, the 
NR spin undergoes elliptical motion with the excursion 
along its e 2 axis (i.e., along the direction of the average 
motion) about several times larger than the oscillations 
along e 3 . The symbols plotted in Fig. 7 reveal that each 
of the elliptic “orbits” corresponds approximately to half 
an orbit of the binary, consistent with the interpretation 
of this motion as nutation. The features exhibited in 
Fig. 8 are similar across all the single-spinning precessing 
cases considered in this work. The small variations in 
spin direction exhibited in Fig. 7 are orders of magnitude 
smaller than parameter estimation capabilities of LIGO, 
e.g. [51], and so we do not expect that these nutation 
features will have a negative impact on GW detectors. 

Let us now apply our nutation analysis to the or¬ 
bital angular momentum directions t. Analogous to 
the spin, we compute averages £ NR and £ PN , and com¬ 
pute the angle between the directions of the averages, 

ZL = Z(>V NR ). This angle—plotted in the top 
panel of Fig. 8 —agrees very well with the difference ZL 
that was computed without orbit-averaging. This indi- 
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FIG. 8. Characterization of nutation effects of the orbital 
angular momentum. Top: angle ZL between the “averaged” 
l in PN and NR for the configuration q5_0.5x (thick red line). 
ZL is shown in thin black line for comparison (cf. Fig. 6). 
The thin blue line shows Z(l, £) between the averaged and 
the filtered signal. Note that it is larger than both ZL and 
ZL. Bottom: the projection of f 1 '* 11 (gray) andt PN (red) onto 
the e .2 — plane described in the text for case q5_0.5x (cf. 
Fig. 8). The system is shown in the interval [—6662, —1556]. 
Both PN and NR show the same behavior, in contrast to the 
behavior of the spin in Fig. 7. The PN-NR matching interval 
is indicated by the horizontal line in the top panel. 


cates that the nutation features of £ agree between NR 
and PN. The top panel of Fig. 9 also plots the angle be¬ 
tween the raw f 1411 and the averaged £ NR , i.e. the open¬ 
ing angle of the nutation oscillations. As is apparent in 
Fig. 8, the angle between and / NR is about 10 times 
larger than the difference between NR and PN (ZL or 
ZL), confirming that nutation features are captured. The 
lower panel of Fig. 8 shows the projection of £ orthogonal 
to the direction of the average £. In contrast to the spins 
shown in Fig. 7, the nutation behavior of £ is in close 
agreement between NR and PN: For both, £ precesses 
in a circle around with identical period, phasing, and 
with almost identical amplitude. We also point out that 
the shape of the nutation features differs between l and 


X±: £ circles twice per orbit around its average £, on an 
almost perfect circle with equal amplitude in the e -2 and 
e-s direction. 

We now extend our precession dynamics analysis to the 
remaining five primary precessing NR simulations listed 
in Table I. The top left panel of Figure 9 shows ZL. The 
difference in the direction of the normal to the orbital 
plane is small; generally ZL < 10 -2 radians, except close 
to merger. Thus it is evident that the trends seen in 
Fig. 4 for ZL hold across all the precessing cases. To 
make this behavior clearer, we parameterize the inspiral 
using the orbital phase instead of time, by plotting the 
angles versus the orbital phase in the NR simulation, as 
shown in the top right panel of Fig. 9. Thus, until a 
few orbits to merger PN represents the precession and 
nutation of the orbital plane well. 

The bottom left panel of Fig. 9 establishes qualitatively 
good agreement for Zyq, with slightly higher values than 
ZL. As already illustrated in Fig. 6, nutation features 
dominate the difference. Averaging away the nutation 
features, we plot the angle Z^i between the smoothed 
spins in the bottom left panel of Fig. 9, where the behav¬ 
ior of Zyi is very similar to that of ZL. This confirms 
that the main disagreement between PN and NR spin dy¬ 
namics comes from nutation features, and suggests that 
the secular precession of the spins is well captured across 
all cases, whereas the nutation of the spins is not. 

All configurations considered so far except 
ql.97_randonr have S ■ £ = 0 at the start of the 
simulations, where S = Si + S 2 is the total spin angular 
momentum of the system. When S ■ £ = 0, several terms 
in PN equations vanish, in particular the spin orbit 
terms in the expansions of the binding energy, the flux 
and the orbital precession frequency, see Eqs. (A14), 
(A15), and (A31) in Appendix A. 

To verify whether S ■ £ = 0 introduces a bias to our 
analysis, we perform our comparison on an additional 
set of 31 binaries with randomly oriented spins. These 
binaries have mass ratio 1 < q < 2, spin magnitudes 
0 < xi ,2 < 0-5, and correspond to cases SXS:BBH:0115 
- SXS:BBH:0146 in the SXS catalog. Fig. 10 plots ZL 
for these additional 31 PN-NR comparisons in gray, with 
ql.97_random highlighted in orange. The disagreement 
between PN and NR is similarly small in all of these 
cases, leading us to conclude that our results are robust 
in this region of the parameter space. 

B. Orbital Phase Comparisons 

Along with the precession quantities described above, 
the orbital phase plays a key role in constructing PN 
waveforms. We use "Pa, a geometrically invariant an¬ 
gle that reduces to the orbital phase difference for non- 
precessing binaries (cf. Sec. IIC) to characterize phasing 
effects. We focus on single spin systems with mass-ratios 
from 1 to 8, where the more massive black hole carries a 
spin of Xi = 0-5, and where the spin is aligned or anti- 
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FIG. 9. Comparison of orbital plane and spin precession for the primary six precessing NR simulations. Top Left: ZL as a 
function of time to merger. Top right: ZL as a function of orbital phase. Bottom left: Z\i as a function of orbital phase. 
Bottom right: Z\i between the averaged spins. All data plotted are averages over 12 matching intervals, cf. Fig. 3, utilizing 
the Taylor T4 PN approximant. The thin horizontal lines in the top left panel show the widest edges of the PN matching 
intervals. 



FIG. 10. ZL for additional 31 precessing configurations 
with arbitrary oriented spins as well as the case ql.97_random. 
Here q € (1,2), xi ,2 < 0.5. For all cases, ZL < 0.5° through¬ 
out most of the inspiral. All data plotted are averages over 
12 matching intervals, cf. Fig. 3. 


aligned with the orbital angular momentum, or where the 
spin is initially tangent to the orbital plane. We match 
all NR simulations to post-Newtonian inspiral dynamics 
as described in Sec. IID, using the 12 matching inter¬ 
vals specified in Eqs. (19) and (20). We then compute 
the phase difference 4 >a at the time at which the NR 
simulation reaches orbital frequency mfl = 0.3. 

The results are presented in Fig. 11, grouped based on 
the initial orientation of the spins: aligned, anti-aligned, 
and in the initial orbital plane. For aligned runs, there 
are clear trends for Taylor T1 and T5 approximants: for 
Tl, differences decrease with increasing mass ratio (at 
least up to q = 8); for T5, differences increase. For Tay¬ 
lor T4, the phase difference 4 >a has a minimum and there 
is an overall increase for higher mass ratios. For anti¬ 
aligned runs, Taylor T5 shows the same trends as for 
the aligned spins. Taylor T4 and Tl behaviors, however, 
have reversed: T4 demonstrates a clear increasing trend 
with mass ratio, whereas Tl passes through a minimum 
with overall increases for higher mass ratios. Our results 
are also qualitatively consistent with the results described 
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FIG. 11. as a function of mass ratio for BBH systems 

with xi — 0-5, and spin direction aligned (top), orthogonal 
(middle), and anti-aligned (bottom) with the orbital angular 
momentum. For clarity, the aligned/anti-aligned data are off¬ 
set by +0.5 and —0.5, respectively, with the thin horizontal 
black lines indicating zero for each set of curves. Plotted is 
averaged over the 12 matching intervals, cf. Fig. 3, and 
for three different Taylor approximants. 

in [52] as we find that for equal mass binaries, the Tay¬ 
lor T4 approximant performs better than the Taylor T1 
approximant (both for aligned and anti-aligned spins). 

For the in-plane precessing runs, we see clear trends 
for all 3 approximants: Taylor T4 and T5 both show in¬ 
creasing differences with increasing mass ratio, and T1 
shows decreasing differences. These trends for precess¬ 
ing binaries are consistent with previous work on non¬ 
spinning binaries [13], which is expected since for S ■ i 
many of the same terms in the binding energy and flux 
vanish as for non-spinning binaries. Overall, we find that 
for different orientations and mass ratios, no one Taylor 
approximant performs better than the rest, as expected 
if the differences between the approximants arise from 
different treatment of higher-order terms. 

C. Convergence with PN order 

So far all comparisons were performed using all avail¬ 
able post-Newtonian information. It is also instructive to 
consider behavior at different PN order, as this reveals 
the convergence properties of the PN series, and allows 
estimates of how accurate higher order PN expressions 
might be. 

The precession frequency m, given in Eq. (A31), is 
a product of series in the frequency parameter x. We 



FIG. 12. Comparison of PN-NR precession dynamics when 
the expansion order of the PN precession equations is varied. 
Shown is the case q3_0.5x. The top panel shows the precession 
of the orbital plane, and the bottom panel of the spin x\ 
(without and with averaging). All data shown are averages 
over 12 matching intervals, cf. Fig. 3. 

multiply out this product, and truncate it at various 
PN orders from leading order (corresponding to 1.5PN) 
through next-to-next-to-leading order (corresponding to 
3.5PN). Similarly, the spin precession frequencies fli .2 in 
Eqs. (2) and (A32) are power series in x. We truncate the 
power series for in the same fashion as the power se¬ 
ries for w, but keep the orbital phase evolution at 3.5PN 
order, where we use the TaylorT4 prescription to imple¬ 
ment the energy flux balance. For different precession- 
truncation orders, we match the PN dynamics to the NR 
simulations with the same techniques and at the same 
matching frequencies as in the preceding sections. 

When applied to the NR simulation q3_0.5x, we obtain 
the results shown in Fig. 12. This figure shows clearly 
that with increasing PN order in the precession equa¬ 
tions, PN precession dynamics tracks the NR simulation 
more and more accurately. When only the leading order 
terms of the precession equations are included (1.5PN 
order), ZL and Z \i are ~ O.lrad; at 3.5PN order this 
difference drops by nearly two orders of magnitude. 

We repeat this comparison for our six main precessing 
cases from Table I. The results are shown in Fig. 13. It 
is evident that for all cases ZL decreases with increasing 
order in the precession equations with almost 2 orders of 
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FIG. 13. Convergence of the PN precession equations for all 
cases in Table I. The evolution was done with the Taylor T4 
approximant at 3.5 PN order. The leading order spin-orbit 
correction is at 1.5 PN order and the spin-squared corrections 
appear at 2 PN order. Each data point is the average ZL over 
PN-NR comparisons performed using 12 matching intervals, 
cf. Fig. 3, with error bars showing the maximal and minimal 
ZL and Zxi of the 12 fits. 


magnitude improvement between leading order and next- 
to-next leading order truncations. A similar trend is seen 
in the convergence of the spin angle Z\i shown in bot¬ 
tom panel of Fig. 13. The angle decreases with PN order 
almost monotonically for all cases except ql.CLtwospins. 
However, this is an artificial consequence of picking a 
particular matching point at mfl = 0.03: as can be seen 
from the bottom panel of Fig. 12 Z\ t shows large oscil¬ 
lations and it is a coincidence that the matching point 
happens to be in a “trough” of xi- 

So far we have varied the PN order of the precession 
equations, while keeping the orbital frequency evolution 
at 3.5PN order. Let us now investigate the opposite case: 
varying the PN order of the orbital frequency and moni¬ 
toring its impact on the orbital phase evolution. We keep 
the PN order of the precession equations at 3.5PN, and 
match PN with different orders of the orbital frequency 
evolution (and TaylorT4 energy-balance prescription) to 
the NR simulations. We then evaluate <I>a (a quantity 
that reduces to the orbital phase difference in cases where 
the latter is unambiguously defined) at the time at which 
the NR simulation reaches the frequency toH = 0.03. We 


examine our six primary precessing runs, and also the 
aligned-spin and anti-aligned spin binaries listed in Ta¬ 
ble I. 

When the spin is initially in the orbital plane, as seen 
in the top panel of Fig. 14, the overall trend is a non¬ 
monotonic error decrease with PN order, with spikes at 
1 and 2.5 PN orders as has been seen previously with 
non-spinning binaries [32]. All of the aligned cases show 
a large improvement at 1.5 PN order, associated with the 
leading order spin-orbit contribution. The phase differ¬ 
ences then spike at 2 and 2.5 PN orders and then decrease 
at 3 PN order. Finally, different cases show different re¬ 
sults at 3.5 PN with some showing decreases differences 
while for others the differences increase. 

For the anti-aligned cases the picture is similar to pre¬ 
cessing cases with a spike at 1 and 2.5 PN orders and 
monotonic improvement thereafter. The main difference 
from precessing cases is the magnitude of the phase dif¬ 
ferences, which is larger by a factor of ~ 5 at 3.5 PN order 
for the anti-aligned cases (see for example ql.5_s0.5x_0). 

These results suggest that convergence of the orbital 
phase evolution depends sensitively on the exact param¬ 
eters of the system under study. Further investigation of 
the parameter space is warranted. 


D. Impact of PN spin truncation 

As mentioned in Sec. IIA 2, post-Newtonian expan¬ 
sions are not fully known to the same orders for spin and 
non-spin terms. Thus, for example, the expression for 
flux T is complete to 3.5 PN order for non-spinning sys¬ 
tems, but spinning systems may involve unknown terms 
at 2.5 PN order; a similar statement holds for dE/dx. 
This means that when the ratio in Eq. (4), J-/(dE/dx), 
is re-expanded as in the T4 approximant, known terms 
will mix with unknown terms. It is not clear, a priori , 
how such terms should be handled when truncating that 
re-expanded series. 

Here we examine the effects of different truncation 
strategies. We focus on the Taylor T4 approximant 
while considering various possible truncations of the re¬ 
expanded form of T/(dE/dx). We denote these possi¬ 
bilities by the orders of (1) the truncation of non-spin 
terms, (2) the truncation of spin-linear terms, and (3) 
the truncation of spin-quadratic terms. Thus, for exam¬ 
ple, in the case where we keep non-spin terms to 3.5 PN 
order, keep spin-linear terms to 2.5 PN order, and keep 
spin-quadratic terms only to 2.0 PN order, we write (3.5, 
2.5, 2.0). We consider the following five possibilities: 

(i) (3.5, 3.5, 3.5) 

(ii) (3.5, 4.0, 4.0) 

(iii) (3.5, 2.5, 2.0) 

(iv) (3.5, 3.5, 2.0) 

(v) (3.5, 4.0, 2.0). 

To increase the impact of the spin-orbit terms, we 
examine aligned and anti-aligned cases from Table I, 
with results presented in Fig. 15. For aligned cases, no 
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FIG. 14. Convergence of the Taylor T4 approximant with PN 
order. Shown are all cases from Table I. Top: all precessing 
cases. Middle: aligned spin cases. Bottom: anti-aligned 
spin cases. Each data point shown is averaged over PN-NR 
comparison with 12 matching intervals, cf. Fig. 3. Error bars 
are omitted for clarity, but would be of similar size to those 
in Fig. 15. 


one choice of spin truncation results in small differences 
across all mass ratios. All choices of spin truncation ex¬ 
cepting (3.5, 4.0, 4.0) have increasing errors with increas¬ 
ing mass ratio. Truncating spin corrections at 2.5 PN or¬ 
der (3.5, 2.5, 2) consistently results in the worst matches. 
On the other hand, we find that, for anti-aligned runs, 
adding higher order terms always improves the match, 
keeping all terms yields the best result, and all choices of 
truncation give errors which are monotonically increasing 
with mass ratio. Overall, anti-aligned cases have larger 
values of $a when compared to cases with same mass 
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FIG. 15. Impact of different choices for spin truncation on 
orbital phase difference 4>a, as a function of mass ratio. The 
lines are labeled by the truncation types, as explained in the 
text. The upper panel shows all cases for which the spins are 
aligned with the orbital angular momentum; the lower panel 
shows the anti-aligned cases. 


ratios. This result is consistent with findings by Nitz 
et al. [14] for comparisons between TaylorT4 and EOB- 
NRvl approximants. 


E. Further numerical considerations 

1. Numerical truncation error 

Still to be addressed is the effect of the resolution of 
NR simulations in the present work. The simulation 
ql_twospins is available at four different resolutions la¬ 
beled Nl, N2, N3 and N4. We match each of these four 
numerical resolutions with the Taylor T4 approximant, 
and plot the resulting phase differences $a in Fig. 16 
as the data with symbols and error bars (recall that the 
error bars are obtained from the 12 different matching 
regions we use, cf. Fig. 3). All four numerical reso¬ 
lutions yield essentially the same $a- We furthermore 
match the three lowest numerical resolutions against the 
highest numerical resolution N4 and compute the phase 
difference 4 >a- The top panel of Figure 16 shows $a 
computed with these 4 different numerical resolutions. 
All the curves lie on top of each other and the differences 
between them are well within the uncertainties due to 
the matching procedure. The bottom panel shows the 
differences in 4 >a between the highest resolution and all 
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FIG. 16. Convergence test with the numerical resolution of 
the NR simulation ql_twospins. Top panel: <I>a with com¬ 
parisons done at different resolutions. All the curves lie within 
uncertainties due to the matching procedure, indicating that 
numerical truncation error does is not important in this com¬ 
parison. The difference between each curve and the highest 
resolution are of order 15% and are within the matching un¬ 
certainties. Bottom panel: ZL with comparisons done at 
all the resolutions. The curves lie within the matching uncer¬ 
tainties. 
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FIG. 17. Gauge change during numerical simulation q5_s0.5x. 
The solid curves represent the recent re-run of q5_0.5x that is 
analyzed in the rest of this paper. The dashed curves repre¬ 
sent an earlier run SXS:BBH:0058 which changes the gauge 
at t — tmerge ~ —3200 M. Top: behavior of the orbital fre¬ 
quency mil in evolution with (dashed curve) and without 
gauge change (solid curve). Bottom: <I>a for all Taylor ap- 
proximants. To avoid matching during the gauge change, the 
matching was done with mfl c = 0.017. 


others. Throughout most of the inspiral, the difference is 
~ 10%. Similar behavior is observed in other cases where 
multiple resolutions of NR simulations are available. We 
therefore conclude that the effects of varying numerical 
resolution do not impact our analysis. 


2. Numerical gauge change 

The simulation SXS:BBH:0058 in the SXS catalog uses 
identical BBH parameters than q5_0.5x, but suffers from 
two deficiencies, exploration of which will provide some 
additional insights. First, the switch from generalized 
harmonic gauge with fixed gauge-source functions [32] to 
dynamical gauge-source functions [41, 42] happens near 
the middle of the inspiral, rather than close to merger 
as for the other simulations considered. This will give us 
an opportunity to investigate the impact of such a gauge 
change, the topic of this subsection. Second, this simu¬ 
lation also used too low resolution in the computation of 


the black hole spin during the inspiral, which we will dis¬ 
cuss in the next subsection. We emphasize that the com¬ 
parisons presented above did not utilize SXS:BBH:0058, 
but rather a re-run with improved technology. We use 
SXS:BBH:0058 in this section to explore the effects of its 
deficiencies. 

While the difference between PN and NR gauges does 
not strongly impact the nature of the matching results, 
a gauge change performed during some of the runs does 
result in unphysical behavior of physical quantities such 
as the orbital frequency. Figure 17 demonstrates this 
for case q5_s0.5x. The old run SXS:BBH:0058 with the 
gauge change exhibits a bump in the orbital frequency 
(top panel), which is not present in the re-run (solid 
curve). When matching both the old and the new run 
to PN, and computing the phase difference <1>a, the old 
run exhibits a nearly discontinuous change in 4 >a (bot¬ 
tom panel, dashed curves) while no such discontinuity is 
apparent in the re-run. 
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FIG. 18. Top: The magnitude of the spin as a function of 
time in the original run (black) and the new run (blue) as 
well as the value computed with the procedure described in 
the text (orange). Middle panel: angles between the spins 
and normals to the orbital plane (thin curves) and their aver¬ 
aged values (bold curves) for the original run and the re-run. 
Lower panel: Z\i and /.l for both the old run and the re¬ 
run (the data of this panel are averaged over 12 matching 
intervals, cf. Fig. 3). To avoid matching during the gauge 
change, the matching was done with mfl c = 0.017. 



q5_0.5x), too low numerical resolution was used for these 
two steps. While the evolution itself is acceptable, the 
extracted spin shows unphysical features. Most impor¬ 
tantly, the reported spin magnitude is not constant, but 
varies by several per cent. Figure 18 shows as example 
Xi from SXS:BBH:0058. For t — f merge < 3200 M oscil¬ 
lations are clearly visible. These oscillations vanish at 
t — fmerge ~ 3200M, coincident with a switch to damped 
harmonic gauge (cf. Sec. HIE 2). Similar oscillations in 
q3_0.5 disappear when the resolution of the spin compu¬ 
tation is manually increased about 1/3 through the inspi¬ 
ral, without changing the evolution gauge. Our new re¬ 
run q5_0.5x (using damped harmonic gauge throughout), 
also reports a clean Xi > cf. Fig. 18. Thus, we conclude 
that the unphysical variations in the spin magnitude are 
only present if both the resolution of the spin computa¬ 
tion is low, and the old gauge conditions of constant H a 
are employed. 

The NR spin magnitude is used to initialize the PN 
spin magnitude, cf. Eq. (14). Therefore, an error in 
the calculation of the NR spin would compromise our 
comparison with PN. For the affected runs, we correct 
the spin reported by the quasi-local spin computation by 
first finding all maxima of the spin-magnitude \ between 
500M and 2000M after the start of the numerical sim¬ 
ulation. We then take the average value of % at those 
maxima as the corrected spin-magnitude of the NR sim¬ 
ulation. Figure 18 shows the case q5_0.5x as well as the 
rerun described in Sec. HIE2. It is evident that this 
procedure produces a spin value which is very close to 
the spin in the rerun where the problematic behavior is 
no longer present. Thus, we adopt it for the three cases 
where an oscillation in the spin magnitude is present. 

The nutation features shown in Fig. 7 are qualitatively 
similar for all our simulations, independent of resolution 
of the spin computation and evolution gauge. When the 
spin is inaccurately measured, the nutation trajectory 
picks up extra modulations, which are small on the scale 
of Fig. 7 and do not alter the qualitative behavior. 

The lower two panels of Fig. 18 quantify the impact of 
inaccurate spin measurement on the precession-dynamics 
comparisons performed in this paper: The middle panel 
shows the differences between the spin directions in the 
original 0058 run and our re-run q5_0.5x. The spin direc¬ 
tions differ by as much as 0.01 radians. However, as the 
lower panel shows, this difference can mostly be absorbed 
by the PN matching, so that Z\i and XL are of similar 
magnitude of about 10 -3 radians. 


3. Problems in quasi-local quantities 

Computation of the quasi-local spin involves the so¬ 
lution of an eigenvalue problem on the apparent hori¬ 
zon followed by an integration over the apparent hori¬ 
zon, cf. [53-55]. In the simulations ql.0_0.5x, ql.5_0.5x 
and q3.0_0.5x and in SXS:BBH:0058 (corresponding to 


IV. DISCUSSION 

We have presented an algorithm for matching PN pre¬ 
cession dynamics to NR simulations which uses con¬ 
strained minimization. Using this algorithm, we perform 
a systematic comparison between PN and NR for process¬ 
ing binary black hole systems. The focus of the compar¬ 
ison is black hole dynamics only, and we defer discussion 
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of waveforms to future work. By employing our matching 
procedure, we find excellent agreement between PN and 
NR for the precession and nutation of the orbital plane. 
The normals to the orbital plane generally lie within 10~ 2 
radians, cf. Fig. 9. Moreover, nutation features on the 
orbital time-scale also agree well between NR and PN, 
cf. Fig. 8. 

For the black hole spin direction, the results are less 
uniform. The NR spin direction x^ R shows nutation fea¬ 
tures that are qualitatively different than the PN nuta¬ 
tion features, cf. Fig. 7. The disagreement in nutation 
dominates the agreement of with Xi* N ; averaging 
away the nutation features substantially improves agree¬ 
ment, cf. Fig. 6. The orbit-averaged spin directions agree 
with PN to the same extent that the £ direction does 
(with and without orbit averaging), cf. Fig. 9. 

Turning to the convergence properties of PN, we have 
performed PN-NR comparisons at different PN order of 
the precession equations. For both orbital angular mo¬ 
mentum t and the spin direction X \, we observe that 
the convergence of the PN results toward NR is fast and 
nearly universally monotonic, cf. Fig. 13. At the highest 
PN orders, the spin results might be dominated by the 
difference in nutation features between PN and NR. 

The good agreement between PN and NR precession 
dynamics are promising news for gravitational wave mod¬ 
eling. Precessing waveform models often rely on the post- 
Newtonian precession equations, e.g. [15, 56]. Our results 
indicate that the PN precession equations are well suited 
to model the precessing frame, thus reducing the prob¬ 
lem of modeling precessing waveforms to the modeling of 
orbital phasing only. 

The accuracy of the PN orbital phase evolution, un¬ 
fortunately, does not improve for precessing systems. 
Rather, orbital phasing errors are comparable between 
non-precessing and precessing configurations, cf. Fig. 14. 
Moreover, depending on mass-ratio and spins, some Tay¬ 
lor approximants match the NR data particularly well, 
whereas others give substantially larger phase differences, 
cf. Fig. 11. This confirms previous work [9, 12, 13, 33, 57] 
that the PN truncation error of the phase evolution is im¬ 
portant for waveform modeling. 

We have also examined the effects of including par¬ 
tially known spin contributions to the evolution of the 
orbital frequency for the Taylor T4 approximant. For 
aligned runs, including such incomplete information usu¬ 
ally improves the match, but the results are still sensitive 
to the mass ratio of the binary (top panel of Fig 15). For 
anti-aligned runs, it appears that incomplete information 
always improves the agreement of the phasing between 
PN and NR (bottom panel of Fig 15). 

In this work we compare gauge-dependent quantities, 
and thus must examine the impact of gauge choices on 
the conclusions listed above. We consider it likely that 
the different nutation features of are determined by 
different gauge choices. We have also seen that different 
NR gauges lead to measurably different evolutions of x, 
£, and the phasing, cf. Fig. 17 and 18. We expect, how¬ 


ever, that our conclusions are fairly robust to the gauge 
ambiguities for two reasons. First, in the matched PN- 
NR comparison, the impact of gauge differences is quite 
small, cf. lowest panel of Fig. 18. Second, the near uni¬ 
versal, monotonic, and quick convergence of the preces¬ 
sion dynamics with precession PN order visible in Fig. 13 
would not be realized if the comparison were dominated 
by gauge effects. Instead, we would expect PN to con¬ 
verge to a solution different from the NR data. 
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Appendix A: Post-Newtonian dynamics 

We consider compact object binary with masses mq 2 
and carrying angular momentum Si <2 . The post- 
Newtonian expressions are most conveniently written us¬ 
ing the following symbols: 


m = 

mi + to 2 , 

(Al) 

v = 

777-1777-2 

(A2) 

9 ’ 

777 z 

5 = 

777 1 — 7772 

(A3) 

777 

s = 

+ S 2 , 

(A4) 

si = 

s-i 

m 2 

(A5) 

Sn = 

S-n 

TO 2 

(A6) 

E = 

777 - 777 -j 

-02- 01 , 

7772 777 1 

(A7) 

<ri = 

E ■£ 

m 2 ’ 

(AS) 







17 




Xs 


Xa 

S 0 


SO 


S • n 


TO 2 

\ (X1+X2), 
7 , (Xi - X 2 ), 


TO - TO - 

- 01 H-D2, 

mi m 2 


A 

TO 2 


(A9) 

(A10) 

(All) 

(A12) 

(A13) 


1. Energy and Flux 

The energy and flux are written as power series 
expansion parameter x = (mil) 2 / 3 : 

E(x) = -\mvx (l + J2k =2 a kX k/2 ) , 

Fix) = fv 2 x 5 (1 + J2k=2 hx k/2 ) . 

For the energy, coefficients are given explicitly by: 


a2 = ~4^l2’ 

14 

a 3 = 26(7 1 + yS/, 

27 19 1 2 , n _ 2 

«4 = -y + y^ - y*' + KXs - Xa - 3[(x« 

,1 


-(Xa-^) 2 ] 


+ (^ - U ){Xs + X 2 - 3[(x« • ^ + (Xa • ^) 2 ]} + *{£ • Xa - 3[(x s • e)(Xa ' *)]}, 


a .5 = 11s/ + 3(5(7/ + z' 


61 10 . 
-TT®/-—0(7/ 

9 3 


a 6 = 


675 

~64~ 


34445 205 


96 


576 


155 


35 


v - v 2 - 1 

96 5184 


,'135 367 29 2 \ ,(21 _ 5 2 , 

a 7 = —-— ^ + — v s/ + (5 — - 39u + -V 07 . 

4 4 12 / V 4 4 


Meanwhile for the flux T\ 


, 1247 35 

62 =-Z4 

336 12 ’ 

g 

63 = 47 r - 4s/ - - 5 ( 7 /, 


44711 9271 65 

64 —--f“ - (7 - 

9072 504 18 


v 2 + 


/ 287 

VllfT 


24 


(Xs-lf 


65 = — 


89 

96 

8191 


Iv 

24 


X s 


287 

1)6 


89 


96 


287 


89 . 


- 12/2 (x a ■ t) + +4u)x a + —6{x s ■ e)(Xa • t) - TX<5(Xs • Xa), 


48 


9 13 

„„„ 7T — “Si — TrOdj + Z2 

672 2 16 


583 272 43 r 

S, + T s< " 


bo = 


6643739519 16 


1712 


69854400 


— 7 T 2 - 

3 105 


856, . 

7® “ Y05 log ( 16a; ) + 


-134543 41 


7776 


67 = 


-167TS/ - —^—(5(7/, 

6 

/476645 6172 2810 


V 6804 189 

16285 214745 


z' 2 I si 


504 1728 


27 
193385 
3024 


(9535 1849 1501 

[_I_7/_7 

v 336 


126 


36 


+ 48" ]v ~ 


6(71 


v 2 I 7 r, 


, , 3485tt 13879tt \ ( 7163t 7 13058377 . r 

= I-— + V ) Si + ( — 7 ^- + V 6(71, 


96 


72 


672 


2016 


48 


94403 2 775 , 

- u" - v 3 

3024 324 


in the 

(A14) 

(A15) 


(A16) 

(A17) 

(A18) 

(A19) 

(A20) 

(A21) 

(A22) 

(A23) 

(A24) 

(A25) 

(A26) 

(A27) 

(A28) 


where 7 e denotes Euler’s constant. 
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2. Precession dynamics 


The evolution of the orbital plane is governed by the frequency vo in Eq. (la), which is defined in terms of two 
auxiliary quantities, 7 = m/r and ai = a ■ £: 


. 3 — v 3 cti + 5si q /9 12 — 65 v 9 

l = x{l + — a: +--- x 3 ' 2 + -—- x 2 + 


30 + 8 ;/ 
9 


si + 2 aiS'j x 5 / 2 


1 + v - 


2203 41n* 


2520 192 


229iy 2 | v 31 
36 + 81 


z 3 + 


60 - 127o - 72v 2 16 - 61v - 16 

12 Sl+ 6 


aiS^j x 7/2 


+x 2 (s 2 - 3(s 0 • /) 2 ) 


ai = : — 7 s n + 3<j n S + x 
m 


2 9v 

s n I - 10 


cr n S I - y - 6 


+ x z 


52v 2 59 k 3\ / 17v 2 73u 3 

—w- + -r- + ^ —— + — + - 


V 6 


3r 4 

(s 0 -£)(s 0 -h). 

m 


(A29) 


(A30) 


Note that we have dropped the pure gauge term — ^ ln(r/rg) from 7. We now have 


vo = 


a; 7 

,3/2' 


(A31) 


The spins obey Eqs. (2) with 


fli — £— 


xi I —35 + 2^ + 3 


10^ — 9 r v 2 5 v 9 

- 5 - 1 - 1 - 

16 24 4 16 


32 


4 _ 48 ~ 


32 


16 32 


m J 


3m 


~(Xi ’ n)h - mix 2 + 3 m|(x 2 • n)h 


(A32) 


The expression for f l 2 is obtained by Xi X 2 , mi O m 2 , 5 o —5 and q +>• 1/q. 


We re-expand the right-hand-side of Eq. (A31), and trun¬ 
cate the expansion for vo and at the same power of x 
beyond the leading order. We refer to the order of the last 
retained terms as the precession PN order. For the ma¬ 
jority of comparisons presented in this paper, we truncate 
at 3.5PN; truncation at lower PN order is only used in 
Sec. Ill C. Note that spin-squared interactions imply the 
lack of circular orbits for generic orientations of the spins. 
We neglect these complications in the present work. 


The norm of a quaternion Q is defined by \Q\ 2 = QQ. 
The inverse of a quaternion is Q~ x = <5/|Q| 2 , which 
means that the inverse of a unit quaternion is simply 
its conjugate. The components of a unit quaternion 
R = ro + r satisfy RR = rl + r- r— 1. Unit quaternions 
are usually referred to as “rotors”. Any rotation can be 
expressed as a rotor, where the rotor acts on a vector v 
according to the transformation law 


Appendix B: Useful quaternion formulas 


v' = Rv R . 


(B3) 


We refer the reader to other sources [47, 60] for general 
introductions to quaternions. Here, we simply give a few 
formulas that are particularly important in this paper. 
First, we introduce some basic notation to be used for 
the four components of a general quaternion Q: 

Q = ( 90 , 91 , 92 , 93 ) = 9o + q. (Bl) 

In this notation, the quaternion conjugate is just Q = 
qo — q, and we note that the product of quaternions is 
given by 

PQ =PoQo -P- q+Poq + qoP + p x q . (B2) 


The form of this expression ensures that v' is a pure 
vector; it has zero scalar part. To see this, we note that a 
quaternion has zero scalar part if and only if its conjugate 
equals its negative, which is true of the right-hand side 
above. We can use this fact, along with p-q= — \{pq + 
qp) and the unit-norm property RR = 1, to see that 
the right-hand side above is indeed an isometry. Finally, 
simple arguments using the cross product can show that 
such a transformation preserves orientation, and since 
the origin is fixed, it is therefore a simple rotation for 
any rotor R. 
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1. Exponential, logarithms, and square roots 

The quaternions are closely analogous to complex num¬ 
bers, except that quaternions do not commute in general. 
One striking example of this analogy is Euler’s formula, 
which generalizes quite directly. If we define the expo¬ 
nential of a quaternion by the usual power series, we get 
for a unit vector u 

exp[0 u] = cos 0 + u sin 9 , (B4) 

which is precisely Euler’s formula with i replaced by u. 
Every rotor R = ro + r can be expressed in this form, so 
it is easy to see that the logarithm of any rotor has zero 
scalar part and is given by 

-» r* If* I 

r := logR = T^T arctan -— . (B5) 

\r | r 0 

It is useful to note that the logarithm of a rotor is parallel 
to the vector part of the rotor. Finding the magnitude 
of r, of course, is just the usual square root of the sum of 
the squares of its components. And the arctan function is 
applied to real values, so we can use standard implemen¬ 
tations of the atan2 function to evaluate it. So we see 
that both the exponential and logarithm of quaternions 
are extremely simple and numerically robust to calculate. 

These formulas can also be used to define general pow¬ 
ers of quaternions. For the purposes of this paper, how¬ 
ever, we only need one particular power of a quaternion: 
the square root. More specifically, given two unit vectors 
u and w, we need the rotor that takes w to u by the 
smallest rotation possible, which is a rotation in their 
common plane. This rotor is given [47] by 

R^u = V-uw = ±—f = , • (B6) 

\/2[l - (uw)o] 

In this expression, u w represents the result of quaternion 
multiplication of the quaternions u and w. ( uw)o repre¬ 
sents the scalar part of this product, so that the square 
root in the denominator is acting on a real number. 
The sign ambiguity is generally irrelevant because of the 
double-sided transformation law for vectors, Eq. (B3). 
However, in certain special applications such as rotor in¬ 
terpolation, the sign must be chosen carefully to be con¬ 
tinuous [47]. 

2. Deriving the frame rotor from l and h 

For both numerical relativity simulations and Post- 
Newtonian evolutions we have data about the positions 
and velocities of the black holes, that can be used to de¬ 
rive the frame rotor Rf, cf. Fig. 2. Given positions of 
the black holes as functions of time, it’s a simple matter 
to calculate their unit separation vector ?z, and then to 
calculate i using 

fll = n x h . (B7) 


Going from t and h to the frame rotor Rf, the idea is to 
first rotate z onto l. This will also rotate x onto some x!. 
We then need to rotate x' onto n, while leaving i in place. 
Of course, the n-x' is orthogonal to £, so we just perform 
a rotation in that plane. This is easily accomplished by 
the following formula: 


Ri = \J-iz , 

(B8a) 

\j ft x Rj^ Ri . 

(B8b) 


Again, the square roots are to be evaluated using 
Eq. (B6). 


3. Comparing frame rotors 

Reference [47] introduced a simple, geometrically in¬ 
variant measure Ra that encodes the difference between 
two precessing systems as a function of time, easily re¬ 
duced to a single real number 4 >a expressing the magni¬ 
tude of that difference. These quantities were mentioned 
in Sec. IIC without much motivation, here we briefly re¬ 
view that motivation. 

In general, we assume that there are two (analytical 
or numerical) descriptions of the same physical system, 
and that we have two corresponding frames Rf a and Rf b • 
To understand the difference between the frames, we can 
simply take the rotation that takes one frame onto the 
other. In this case, the rotor taking frame A onto frame 
B is 

Ra '■= RfB RfA- (B9) 

Rotors compose by left multiplication, so it is not hard to 
see that this does indeed take RfA onto RfB because the 
inverse of RfA is just its conjugate, so Ra RfA = Rib- 

A particularly nice feature of Ra is that it is com¬ 
pletely independent of the inertial basis frame (x, y, z) 
with respect to which we define the moving frames. That 
is, if we have another basis frame {x ', y ', z r ). there is some 
Rs such that x' = RgxRg, etc. The frame rotors would 
transform as RfA H > RfA = Ria Rs , in which case we 
obtain 

Rf B R!f a = RfB Rs Rs RfA = RfB RfA- (BIO) 
That is, Ra is invariant. 

Now, we seek a relevant measure of the magnitude of 
the rotation Ra- We know that it may be written as a 
rotation through an angle 4> about an axis v. Clearly, (f> 
is the measure we seek. The rotor corresponding to such 
a rotation is given by R = exp[</>t)/2]. Thus, to find the 
angle, we just use the logarithm: </> = 2|logR|, where the 
norm is the usual vector norm. Again, the formula for 
the logarithm of a rotor is a simple combination of stan¬ 
dard trigonometric functions applied to real numbers, as 
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shown above. Using this interpretation with our differ¬ 
ence rotor, we see that the appropriate definition is 

$A := 2 |log [i?fB Ria] | ■ (Bll) 

There is information contained in the direction of the 
logarithm. For example, the component along £ is re¬ 
lated to the difference in orbital phase for non-precessing 
systems, while the component orthogonal to £ is related 
to the direction and magnitude of the difference in £ it¬ 
self. For the sake of simplicity, however, we focus on the 
magnitude of the logarithm, as given above. 

4. Inadequacy of 4 >a — &b for comparisons of 
precessing systems 

We claim that it is impossible—when analyzing pre¬ 
cessing systems—to compare two rotations Ra and Rb 
in a non-degenerate and geometrically invariant way by 
defining some phases <J>a and <f>£ for them separately, 
and then comparing them as Ta — 4>b. Here, “non¬ 
degenerate” means that the phase difference is zero if 
and only if Ra and Rb represent the same rotation, and 
“geometrically invariant” means that the result is not af¬ 
fected by an overall rotation of the basis used to define 
Ra and Rb- In this section, we prove that statement. 

We begin by defining a function $ such that §(Ra) = 
d >a and $(i?s) = The domain of this function is 

a rotation group, which could be the one-dimensional 
group U(l) for non-precessing systems, but must be the 
full three-dimensional group 8 SU(2) for general precess¬ 
ing systems. The range of $ is the usual range of phases, 
the additive group of real numbers modulo 27r. It will be 
useful to note that this is isomorphic to U(l). Finally, 
non-degeneracy is the condition that <J>a — =0 [or 

equivalently &(Ra) = if and only if Ra = ±Rb- 

The condition of geometric invariance can be written 
as a condition on $ itself. If, for example, we measure 
everything with respect to some basis ( x,y,z ), and then 
measure again with respect to some other basis (x', 
we should get the same answer. Now, there is some rotor 
Rs that takes the first basis into the second. If Ra is 
defined with respect to the first basis, then the equiva¬ 
lent quantity will be Ra Rs with respect to the second. 
Geometric invariance is then the statement 

${Ra Rs) - *{Rb Rs) = ${Ra) - ${Rb), (B12) 

for any choice of Rs in SU(2). We will show that there 
is no such $ because the rotation group SU(2) is not 
isomorphic to U(l). 


Even though it is a double cover of the physical rotation group 
SO(3), we use SU(2) here for consistency of notation, because it 
is the group of unit quaternions. The proof would actually be 
slightly simpler for SO(3); we would have &(Ra) = &(Rb), if 


Since Eq. (B12) is true for any rotor Rg , we can choose 
Rs = Rb 1 , and find that 

<S>(Ra Rb 1 ) - $( 1 ) = ${Ra) - $(Rb). (B13) 

Now, we define another function d>'(i?) = < f > (i?) — <E>(1). 
The last equation becomes 

$'(R a Rb 1 ) = &(Ra) ~ &(Rb)- (B14) 

In exactly the same way, we can see that 

&(Rb Ra 1 ) = &(Rb) - &(Ra) = -&(Ra Rb 1 )- 

(B15) 

This must be true for all values of Ra and Rb, so we 
have shown that 

^'(R- 1 ) = -$'(R), (B16) 

for arbitrary R. Therefore, we can also see from 
Eq. (B14) that 

&{RiR2) = + $'(i? 2 ), (B17) 

for arbitrary R\ and Ro- This is precisely the statement 
that $' is a homomorphism [from SU(2) to the additive 
group of real numbers modulo 27r]. 

However, now we can impose the condition that — 
= 0 if and only if Ra = ±Rb- Using the properties 
of homomorphism, it is clear that this is equivalent to 
the statement that the set of all elements that map to 
0 under d’ 7 is ker<F = {—1,1}. Then, the First Group 
Isomorphism Theorem [61] says that the image of $' is 
isomorphic to SU(2) modulo this kernel, which of course 
is just SO(3). But the image of 4>' is (possibly a subgroup 
of) the group U(l), which is obviously not isomorphic to 
SO(3). 9 Therefore, it is impossible to construct a func¬ 
tion fulfilling our requirements for precessing systems. 

It is, however, interesting to note that if our rotation 
group were not SU(2), but the one-dimensional rotation 
group U(l), there would be no contradiction. This is 
how it is possible to construct a useful measure of the 
form — 4 >b for non-precessing systems, because the 
rotations can be restricted to rotations about the orbital 
axis. On the other hand, for precessing systems, the mea¬ 
sure described in Secs. IIC and B 3 is able to satisfy 
both key features of a useful measure (non-degeneracy 
and geometric invariance) because it simply does not 
attempt to define a homomorphism from the rotation 
group; rather, it defines a (non-homomorphic but invari¬ 
ant and non-degenerate) function from two copies of the 
rotation group onto phases, SU(2) x SU(2) —> U(l). 


and only if Ra = Rb-> an d ker^' = {1}. 

9 This statement will not be controversial, but for form’s sake we 
can prove it simply by noting that U(l) is commutative, whereas 
we can find elements of SO(3) that do not commute. 
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